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LES, DNS and RANS for the Analysis of High-Speed 

Turbulent Reacting Flows 


V. Adumitroaie, P.J. Colucci, D.B. Taulbee and P. Givi 
Department of Mechanical and Aerospace Engineering 
State University of New York at Buffalo 
Buffalo, New York 14260-4400 


Abstract 

The purpose of this research is to continue our efforts in advancing the state of knowledge in 
large eddy simulation (LES), direct numerical simulation (DNS) and Reynolds averaged Navier 
Stokes (RANS) methods for the computational analysis of high-speed reacting turbulent flows. 
In the second phase of this work, covering the period: August 1, 1994 - August 1, 1995, we have 
focused our efforts on two programs: (1) Developments of explicit algebraic moment closures for 
statistical descriptions of compressible reacting flows, (2) Development of Monte Carlo numerical 
methods for LES of chemically reacting flows. This report provides a complete description of 
our efforts during this past year as supported by the NASA Langley Research Center under 
Grant NAG-1-1122. 


Technical Monitor: 


Dr. J. Philip Drummond (Hypersonic Propulsion Branch, NASA LaRC, Mail Stop 197, Tel: 804- 
864-2298) is the Technical Monitor of this Grant. 


1 



1 Introduction 


We have just completed our Year 2 of the Phase II activities on this NASA LaRC sponsored project. 
The total time allotted for this phase is three years; this phase was followed at the conclusion of 
Phase I activities (also for three years). Thus, in total we have completed five years of NASA LaRC 
supported research and one more year is remaining. Within the past five years, we have considered 
many issues of interest to the NASA LaRC in improving the state of affairs in DNS, LES and RANS 
of high speed turbulent reacting flows. Our previous four annual reports provide detail information 
on our past achievements. This report provides a complete description of our activities in Year 5. 

Our efforts within the past year have been primarily concentrated on two main tasks: (1) Develop- 
ment of algebraic moment closures for statistical description of (highly) compressible flows, and (2) 
Development of an efficient Monte Carlo computational procedure for LES of chemically reactive 
flows. The efforts in (1) are in continuation of our previous work 1 (discussed in our Year 4 annua] 
report), and the work pertaining to (2) is in continuation of our previous work 2,3 (discussed in our 
Year 3 annual report). In addition, we have devoted a portion of our efforts to make use of the 
models in (1) for the purpose of LES. At this point, our achievements are not yet a level suitable 
for documentation. Our achievements on each of the two constituents of the work in Year 5 are 
described in the next two sections. 


2 Algebraic Turbulence Closures for High Speed Turbulent Flows 

2.1 Introduction 

For the incompressible regime the literature on computational prediction of nonreactive turbulent 
transport is abundant with schemes based on single-point statistical closures for moments up to 
the “second-order”. 4-8 Referred to as Reynolds stress models (RSM), these schemes are based on 
transport equations for the second order velocity correlations and lead to determination of “non- 
isotropic eddy-diffusivities.” This methodology is more advantageous than the more conventional 
models based on the Boussinesq approximations with isotropic eddy diffusivities (such as the k — e 
type closures 4,9 ). However, the need for solving additional transport equations for the higher order 
moments could potentially make RSM less attractive, especially for practical applications. For 
example, it has been recently demonstrated 10 that the computational requirement associated with 
RSM for predictions of three-dimensional engineering flows is more than 80% higher than that 
required to implement the k — t model. The increase is naturally higher for second-order mod- 
eling of chemically reacting flows owing to the additional length and time scales w’hich have to 
be considered. 11-17 A remedy to overcome the high computational cost associated with RSM is to 
utilize “algebraic” closures. 18-25 Such closures are either derived directly from the RSM transport 
equations, or other types of representations 26-29 that lead to anisotropic eddy diffusivities. One 
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of the original contributions in the development of algebraic Reynolds stress models (ARSM) is 
due to Rodi . 20 In this work, all the stresses are determined from a set of “implicit” equations 
which must be solved in an iterative manner. Pope 18 offers an improvement of the procedure by 
providing an “explicit” solution for the Reynolds stresses. This solution is generated via the use 
of the Cayley-Hamilton theorem, but is only applicable for predictions of two-dimensional (mean) 
flows. The extension of Pope’s formulation has been recently done by Taulbee 19 and Gatski and 
Speziale . 24 In these efforts, the Cayley-Hamilton theorem and the “symbolic” matrix manipulation 
techniques are used to generate explicit algebraic Reynolds stress models which are valid in both 
two- and three-dimensional flows. 

In recent years the fundamental research on compressible turbulent flows has experienced a period of 
impetus owing to an increasing involvement of the propulsion community in the design high-speed- 
high-altitude ramjet engines. Although new experimental and numerical information is continuously 
accumulating over the years (for reviews see Lele , 30 Gutmark et al. 31 ) the theory of compressible 
turbulence has not reached maturity yet. Several important aspects have been recognized about the 
nature of the turbulent state of a compressible medium and progress has been made in advancing 
the modeling of simple physical flow phenomena, but the inclusion of compressibility effects and of 
variable inertia effects in the turbulence models is an issue still under investigation, especially for the 
second-order moment closures. Using dimensional analysis in physical space 32 or in Fourier space , 33 
asymptotic analysis , 34 rapid-distortion theory , 35 singular perturbation method 36 inside acoustic 
theory previous contributions have exploited the decomposition concept of the compressible field 
to generate models for terms specific to high-speed flows, i.e. pressure dilatation and dilatational 
dissipation which have been perceived to contribute to the reduced growth rate of the compressible 
mixing layer. These models have been applied in many instances as compressibility corrections in 
conjuncture with the standard k-( model 37 or with a generalized form of the k-e model 33 or with 
the actualized incompressible Reynolds stress turbulence model . 38,39 By contrast true compressible 
second-order modeling attempts 35,40 are very few. 

The specific objective is to provide explicit algebraic relations for the Reynolds stress and for 
the “turbulent flux” of scalar variables in the high-speed regime. Both non-reacting and reacting 
flows with heat release are considered. In the latter, a second-order irreversible chemical reaction in 
considered in turbulent flows with initially segregated reactants. The closures explicitly account for 
the influence of the the turbulent Mach number and Damkohlernumber (only in the scalar model) 
and density gradient, pressure gradient and mean dilatation effects are included in the closures. 
Similar to previous contributions , 18 - 4119 " 25 the starting equations are the differential equations for 
the second order moments. Linear closures for the pressure-strain and the pressure-scalar gradient 
correlations are proposed and simple models for the averaged Favre scalar fluctuations are derived 
and embedded in the final explicit algebraic models. 
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2.2 Governing Equations 


In the statistical approach to the turbulence problem the instantaneous equations are used to 
obtain the governing equations for the mean variables. Denote by overline ensemble average and 
by brackets density weighted (Favre) ensemble averaging: 



Accordingly, we have the following decomposition rules: 


X = X + A", X' = 0 
A' = (X) + A”, (X”) = 0, JC = X - (X) 


The governing equations are written in normalized form (with respect to reference values: 6 W - 
vorticity thickness for length, p Q c , ««>, Too, Poo) for a compressible, reacting with heat release 
{A + tB —* (r + 1 )P + heat ) mean turbulent flow. 

Continuity: 

? + W =0 . (i) 

dt dx 3 v ’ 

Conservation of momentum: 


dpjui) &p(uj)(uj) _ dp da Jt (u) 

dt dxj ~ dxj dxi dxj 

i,j =1,2,3 (2) 


where the stress tensor is a Jt (u) = 2p[Sij{u)-\S pp (u)f> tJ ]/Re = 2 pS^(u)/Re. The present notation 
for strain rate S tJ (u) = (|f i + |^ i )/2 and for all the other linear differential operator is more suitable 
for the compressible problem where the two type of averages (which have different properties 
with respect to the linear differential operators) are naturally encountered. Hereinafter the star 
exponent will indicate the traceless tensor (deviatoric part) correspondent to the unstarred tensor. 
The mean viscosity (p.) follows a Maxwell-Rayleigh variation law with the mean temperature, i.e. 
(p)/prej = ((T)/Tre/) n , n = 0.76. The fluctuations of the viscosity are neglected so that its 
correlations with other variables in the flow are zero. Within the present notation the averaged 
stress CTji(u) is equal to a Jt ((u)) + a ji(u ). 

Let e t = (T - ~fCeY P + u 3 u } /2) l{^{^ - 1 ) A/ 2 ) . Then the total energy equation is: 


&p{e t ) dp{e t ){u 3 ) _ dq 3 (T,Y P ) dp{u)e t ) d 

~dT + dxi - d^r a*/ J P j) 


( 3 ) 
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where the averaged heat flux is 




Vj(T,Y P ) = Qj((T), ( Y P » + qj (T\Y' P ), 

Wj^PM + ^jpPiujT"), 

u^ji(u) = (Uijffji ((«)) + + wTaji(u”) 


and 


, - 1 «m"\ Ce - - w v , 

««i e . 1 = 7 ( 7 -l)M^ ( “j r 1 • (7- + >W + 2~ ' 


Conservation of species: 


mYo) i WM - \ d ( ^ i ir 

dxj dxj ySciie dx-j J 

a = j4, P. 


H) 


where u> a represents the rate of chemical reaction (u>a — ~^B — — t+y^p): 


u! a — -Da ex p 




(5) 


and its mean approximated as: 


di Q ~ - Da ex p 




P 2 (Y a Y B ), 


( 6 ) 


And in the assumption of a perfect gas mixture, the equation of state: 


V = 


1 

7 M - 2 


P<r> 


(7) 


Here />, u;, p, e ( , T, Y Q , Pe, Pr, Af, Ce, Pe, Sc, Ze and Da denote the fluid density, the I'-th 
component of the velocity vector, the pressure, total energy, temperature, mass fraction of species o, 
the Reynolds number, the Prandtl number, the Mach number, the heat release parameter, the Lewis 
number, the Schmidt number, the Zel’dovich number and the Damkohlernumber, respectively. 

The closure problem consists in providing models or closed transport equations for the second order 
moments that appear in the equations for the mean variables. In general, the models are based on 
length scales and time scales constructed from second order quantities obtained from the following 
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transport equations. The equation for the kinetic energy of the turbulence (k) = («, «, )/ 2 reads: 

d 


dp(k) dp(k ){uj) 


dt 


dxj 


dx 3 


\p{u'jk) + p'u] - uloji(u) 


„ » K d(rii) , du” — dp , —da Jt ((u)) _ 


dij 


dij' 


(8) 


The turbulent dissipation p e = = p(1 s + cj) = 2pSlij(u ) + |//5^ p (u ) is split 

into solenoidal and dilatational parts, 42 ' 32 ’ 34 with (l,j(u) = ±(§^ - |^-) as the rotation rate. The 
solenoidal part is computed from the equation: 


dp dp (j(uj) d 

dt dxj dxj 


p( u 'j e s) - 


p de s 
Re dxj 




CtlP ( dxj 


f* — S 

L/ « 2 r j i.v ? 


I 


( 9 ) 


with C (1 = 1.44 and C (2 = 1.92. For the dilatational part there are available several models. 32-34 - 36 
The model of Taulbee and VanOsdol 33 for the dilatational terms combined requires additional 
transport equations and our intention this work is to keep the number of equations at minimum. 
Other options are the model proposed by Sarkar el a/. 34 


U = <s M ? 


( 10 ) 


and the model of Ristorcelli 


:36 


Q = { [/f + 6/*/|] + (|) ' 5 r 2 [3a 2 + 5a; 2 ] 

) 2 [l3a 2 +15u; 2 ]r 2 a 2 /[ } ^< 


*r + 

b 3 + V 15 


( 11 ) 


The parameters are 7j = 0.3, /j = 13.768, 7J — 2.623, 7[ — 1.392, 7^ — 3, a — 1 I 4. Also, Alt 
denotes the turbulent Mach number, that is M ( 2 = §(fc)/c 2 and R t the turbulent Reynolds number 
R t = p(k) 2 /(9e 3 p)Re. The local speed of sound is given by c 2 = T/M 2 . 

The corresponding second-order quantities necessary for the temperature calculati ons are the tem- 
perature variance (T" 2 ) and the thermal turbulent dissipation, pe$ = • F° r lh e former, 

starting from the temperature equation with c v constant, we have: 


dp(T" 2 ) dp(T" 2 )( Uj ) dp(u]T ” 2 ) d ( jp dr 2 

^ — Ci I Q.. n-D. 


d(T) 


dt 


dxj 


OX: 




2-)p dT 


PrRe dx 


2(7-1 )p[(T)(T 


yjl + ± (22 >lt 8 JD) _ ivLvr&p + 27Ce -y, 

•j dxj dxj yPrRe dxj J PrRe dxj dxj 


, du 'j 

dx. 


+{T -^_M + / T yju 


dxi 

+ 27 ( 7-1 )M 2 T"a ij (u)^. 


(12) 
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The transport equation for the thermal turbulent dissipation is given by: 


&p €,> , &pe^(uj) d x , M ^ 

IT + - fe— = ~Wi + s7FeW, 


d (_, » v , H dc#\ ^ _ (s ! ^^ K d{T) 
dxj + ScRedxj) CyiP (k) U i T ) dxj 


CyiPjTX^i u ]) 


d(u,) 


S-1 — S~1 — ^ 

(fcj x “• “> ' dxj 2\~^y*P m • 


<*> 

’ {T" 2 ) 


W (13) 

The constants take the vaJues C yi = 2.0, C y3 = 2.0, C Vi = C c 2 — 1 = 0.92, C y2 = 0.5. 

Treatment of the scalar variable requires the solution of addition al transport equations for the 
reactants’ covariance (Y"Yp) and dissipations pe Q p = For the covariances we have: 

mY'cYj) = _3_ 

dxj 5xj l Sci?e 


. - . X gOi 

1 dx^ dxj 

- v .. w) v -»^> _ 2/r by^y; 

dx, dx , ScRedx , dx. 


+w Q y; + ^y;. (h) 

Source terms in the expanded form read (no summation on greek indexes in all subsequent 
tions) 


equa- 


oJ Q Yp + upY" = -Da exp 


-Ze 


1 1 


,(F) ryy J 


? 2 K ( 07 > + wmyb) + nr: y' b ) + 

( v;r B )){Y A ) + (>;»«> + (i'a'y;yg)]. («> 

Similarly, the scalar dissipations are obtain from (hereinafter c a — ) • 


&pe o0 dp € a p(uj) _ d_ . P 

at + axj azj t P J ap ScRc 9x j ) 

/ a ... i n i. . t \ 


„ 1 (, .„. v 9(y„) , 8{«>) 

- c «' , (t)2 p’P-fa” + - c » p Tk) { ' J> ^ 


’ (k) x “* “ J/ dx, 

n — — £f _£ap c 

” (y;y;> ^ <*> ° fl ' 

In this equation, the chemical source term is of the form: 

5c,/? = — Da exp 


(16) 


Ze (<T) 7y 


J P 2 [(C/la + €A0)(Yb) + ( € B0 + €«B)(V^)] • 


(17) 


To close the transport equations for the second order quantities all the third-order transport terms 
are described by the gradient diffusion hypothesis. Denoting by a"b" any of the second-order 
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correlations, we have: 


(18) 


r -(k) . « d(ab "> 

p{ u i a t> ) = g » 

where C 8 is taken to be equal to 0.22 for all non-gradient correlations ( a b” = k or ab" = V^ 2 ), 
whereas for the turbulent dissipations ( a b = c s or a b = e a ), C 5 = 0.18. Also, the molecular 
transport terms are neglected under the assumption of high Reynolds-Peclet numbers flow. 


2.3 Models Development 


2.3.1 ARSM 


An improved explicit ARSM has been derived by Taulbee 19 from the modeled transport equation 
for the Reynolds stresses. This model is based on the general linear pressure-strain closure given by 
Launder et a/. 43 The improvement is due to an extended range of validity; the model is valid in both 
small and large mean strain fields and time scales of turbulence. A similar line of reasoning is made 
to obtain an algebraic closure for the unclosed correlations in compressible regime. The transport 
equations governing these correlations are transformed into algebraic expressions by making two 
assumptions: (1) Existence of a “near-asymptotic” state, and (2) the difference in the transport 
terms is negligible, in other words we look for the fixed point solution in the structural equilibrium 
limit. The starting equations for the Reynolds stress equation are described by: 


dp{ u ”i M j) d_ 

dt dx k dx k 


[/)(«’’«” u k ) + p'u"6 t k + p'u" Sjk - u’Jatjfu”) - u”CTfcj(u”)] 


-p(u’u k ) 


djuj) 

dx k 



— cfp 
U > 0x t 


—da ki ((u)) 

+u ’-d^r 


. 7 r da k j{(u)) 
‘ dx k 


„ On „ du] 

°‘ k(u ] dt k + ° Uu ) aT k 


(19) 


Hereafter, a designates the anisotropic stress tensor, a,j = [(u t - tij)/(fc) — 26 tJ /3], the Kronecker 
symbol is 6 = [<5 tJ ] = 1 for i = j = 1,3 and 0 otherwise, r = ( k)/e s is the local turbulence time 
scale, o = (5* J ((u))5* t ((u))) 1/2 and w = (ft,j({w))ftj;({u))) 1/2 are tensor invariants. 

The Reynolds stress equation is rewritten in terms of a/(r<r): 


_Da, J /(ro) 

TOp 


(k) 


m r dT kl - 
[dx k dx k . 


Dt 

1 

" W) 


1 

W) 


dTjjk _ ( u "> u ])dT k 
dx k ( k ) dx k 


P* } + n* + M*j + V tJ 


* P<i I- 




\c, 


1 2 


2 + (2 -c., )= L + I^ + 2 m + v+^-v,. 

p e s ° Dt p e. 


( 20 ) 
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where 


Pii = Pii - o Ptii = -P 


, n » d(llj) , , - » K 0\Ui} X, - n U\Ulf | 

airM 


#<«.} 


d<M/> 


( 21 ) 


is the production of Reynolds stress, 


n*, = n tJ - -p'dSij = p' 


(duj_ 

\dx-j dxi ) 


2 du, 

3 P d^ k 


v t j 


is the pressure-strain correlation, 

... ,, 2 -»«»«((«» — ^(W), 

Vo = v„ - -V^ = «, +«, . - 3 », 9 K, 


3 ' ’ 3 dx k 

is the mass flux/ viscous diffusion term, 


M*j = Mi 


3 MSij = 


— dp —dp 2— dp 

Ui dT : +u >frT,- 3 Uk W k ' 3 


(22) 


(23) 


(24) 


is the mass flux/pressure gradient term (also called enthalpic production by exchange with enthalpic 

42 \ 

energy ) 


P «ii = P\ e V ~ o (6 


, „ . du” . „du” 2 „ du", 

= °* {u fe + ° k ' (u _ s' 7 ' 1 *" W” 


(25) 


the anisotropy of the dissipation. The dyads involving the mass flux vector can be added to a single 
tensor. In free shear flows the viscous diffusion part is negligible owing to the high Reynolds numbers 
characteristic to these flows. Nevertheless, the present analysis can accommodate the discarded 
term when necessary, such as near- wall flows. Hence, in subsequent equations the tensor M tJ can be 
used to include the viscous effects. Using a rationale similar to the incompressible situation, 43,6 the 
pressure strain- correlation model can be written as (in this attempt the supplementary compressible 
terms that appear in the Poisson equation have been neglected as being of second order, a rough 
approximation, but valid in the low Mach number regime): 


n tJ - p c* = A , j + 2p(l p i g j + I pjgi )(S pq ((u)) + ft p ,«tt)) (26) 

Our goal being to obtain an explicit algebraic Reynolds stress model, the integrals Ipiqj as well as 
the tensor A, 3 will have to be expressed as linear functions of the anisotropy of the Reynolds stress 
tensor so that the final equation is solvable by exact analytic methods. Therefore 


■'PWJ 


(*) 


Aij — Ci p cnjj T App 

= CtlSgiSpj + a 2 (^pq^tj + bqjbpi) + 0,\b p] a qi + 
a 2{ftpq a ij + fipiQq] T J a pq + ^jq a pi) T a 4&qi a p] 


(27) 
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This above form for l piqj satisfies already symmetry constraints. To determine the coefficients two 
more constraints are applied: the normalization condition which translates into l ppqj - 
and a matching condition such that from the trace of 11* is obtained an existing model for the 
pressure dilatation (replacing the customary incompressible constraint which is recovered from the 
above condition in the limit of zero Mach number). Some of the existing proposals to model this 
term necessitate transport equations such as density variance 33 or pressure variance. 35 Two recent 
pressure dilatation models 44 ' 36 do not require separate equations, the model of Sarkar: 44 


In 

2 


pp 


— —p'd = — YiM?[— (-- — )p(/:)5'pp + Y 7 P f s] 

- pa- 3xi Xi 


(28) 


where X i = 0.15, X 2 = 0.2, X 3 = 0.2 and the model of RistorceUi 36 

Td = -xM?[-\l>(k)S pp + P - p e + T k - 7(7 - l)(Pr + P 7 + Tt)} 

, £(3o 2 - 5d 2 ) 

p(k)M t \ 


(29) 


where 

J r 

2l vd f 1 pd _____ 

X = 1 + 2 ~Y X = l + HpdMi + ilrdMMi-l)' 

I P d = | 7 i + F vd [ 2cr2 “ ' I P d= 30 (3) ° 7l 

Using the latter model, the matching constraint on the I piq] produces 

Ipiq ,i/(&) = ly L ( g^P9 + a p?)- 


(30) 


The same condition enforces also that 

A pp = X^tiP 7 + 7 


, , D(3ct 2 — 5ct 2 ) 

1 ){Pt + P { )1 — pfi)M t X 


(transport terms have been neglected based on the local homogeneity assumption). In this manner 
a linear pressure-strain model is obtained. It is known that with linear forms it is impossible 
to satisfy realizability conditions requiring that the eigenvalues of the Reynolds stress tensor be 
positive. To overcome this deficiency we employ a method suggested by Schumann 45 and detailed 
by Shih and Shabbir. 46 If F = 1 + 27/77/8 + 9/7/4 is a parameter involving the second invariant 

// = -1 and third invariant 777 = of the Reyn ° ,dS Stf€SS anis ° tr ° py tenS ° r ’ 

then the following asymptotic behavior for the pressure strain-model ensures that realizability is 

satisfied: 


A 


ee 



as F 


0 
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d(up ) 7 

~dT lpeq£ 


0 as F — ► 0 


(31) 


where the index e indicates that the relations are written in the principal axes of (*<«,■)• To enforce 
this kind of decay additional parameters are introduced in the pressure strain-model which reads 

in final form: 

n *j - p e *a = - c '~p ?a 'j- 4rF “ r + pW [[I + ^ xM ‘ 2 ] 


[1 - C 3 + xM?) a ip S pj ((u )) + s* p ((u))a pj - -5^((u))a p ^ t jj 
[l-C 4 - xM?}[a ip ftpj({u)) - Q,p((u))a P j] + ^ X^ 2 •S'ppC < ta ) ) a ii] B r F 0r 


(32) 


with C 3 = (5-9 C 2 )/U and C 4 = (1 + 7C 2 )/11. The value for the constant C 2 will be the same as in 
the incompressible model to preserve consistency in the zero Mach number limit, that is C 2 - 0.45. 
The parameters are a, = 0.1, fir = 0.5, A r = min(F— , 0.1—) and = min(f*, 0.1' '). The 
mass flux is usually modeled by gradient-transport hypothesis 38 or solving its transport equations.^ 
A compromise between economy and accuracy is obtained using a model proposed by Ristorcelh. 


P u . = 


= T„ 


dM . 7 d{ui)d{u k ) 
v 0 bij + v x t u — r v 2 t u 


dx t 


dxk dxj 


i> .t &p 

( u : u v)q x 


(33) 


where r u = M t r/[ 1 + <) - !)]• Here t/ 0 , and ^ 2 are the coefficients obtained from 

the inversion of the matrix G tJ = 6 i} + r u d 4^: v 0 = -(1 + 7 G + /7 g)^, ^ - J 1 + <■)«*» 
v 2 = (1 + / g + // G + 777 g)" 1 » tlie Roman llumbers ^Presenting the invariants of G. For simplicity 

we will use only the lowest order contribution from this model to obtain: 

M*j = - = |^Mo(*)(^ij + n i' - \ n rP 6 n) = -^ T ^o( k ) R >: ^ 


where _ dp &p_ 

dxi dxj 

is the mean density gradient-mean pressure gradient dyad. If the entire model is to be used in the 
following equations R* should be replaced with M* and set b 2 = -1 /(/>(*»• 

The final equation is obtained introducing the above expressions for the pressure-strain correlation 
model, production and mass flux terms into the equation for a*. To get the fixed point solution we 
set to zero the Lagrangian derivative as well as the difference m the transport terms. us res 
in the following linear tensorial equation written in matrix form (the braces signify the trace o e 

matrix inside the braces): 


a = -gr 


bi S* + b 2 R* + h faS* + S*a - ^{aS*}tf) - M aff Ha)j 


(35) 
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with h = i-B r F^(i + ^ X M?),b 2 = j,lT u v 0 ,b 3 = l-B r F0r(\-C 3 + X M?),b 4 = l-B r F^{\- 
C 4 ~ XM?), and 


9 = 


r 6 P T Dcr 

{ArF^C,- + C (2 - 2 + (2 - C £l )— + - — + 


P <-> 


^(1 - xM 2 B t F 0r )S pp ((u)) + 2 — + ^ 

O P 


l-l 


(36) 


The task of solving the equation (35) is formidable. From the above equation we see that the 
anisotropy of the Reynolds stress tensor is dependent on three primary second order tensors , two 
symmetric and one skew-symmetric, a,y = a t j(S*, fi, R*). The solution can be expressed as a finite 
3-D tensor polynomial, 


a = £C A T A (37) 

A 

that is a linear combination of all the linear independent tensor products formed from the three 
primary tensors. The coefficients of this polynomial are function of the set of independent invariants 
which form the integrity basis for this problem. In this case the dimension of the minimal tensor 
base is A = 41 (cf. Spencer 48 ) and although not all of the tensor products will appear in the final 
result, there is little hope that the model will be of practical use at this time. A way to circumvent 
this difficulty is to make a simplifying approximation regarding the b 3 term. For most practical 
free-shear flow applications F > 0.1, therefore in the low turbulent Mach number domain b 3 ~ C 3 . 
It has been argued 19,25 that for the range of values used for the constant C 2 the inequality C 3 <C C 4 
holds and therefore the term multiplied by C 3 will have small effect on the solution. Thus using 
the superposition principle a = a~' + a 1 ^ where a^ stands for the solution dependent on S*, 

a s = -gr [fqS* - 6 4 (a s fi - tta s )] (38) 

and denoting the solution dependent on R*, 

a R = - gT [l> 2 R* - b 4 {a R n - Oa R )] . (39) 

Applying the results of Taulbee 19 we have 

a s = — 2 C m tS* - 4a 2 r 2 (S*ft - CIS*) - 8a 3 r 3 (fi 2 S* + S*ft 2 - ^{S*fl 2 }£) 

-16a 4 r 4 (fiS*fi 2 - fl 2 S*fl) - 32o 5 r 5 {S*fi 2 }(n 2 - ^{fi 2 }£) (40) 

where = b 4 g(l — |/ioCT 2 )/ii, a 2 = \b\b 4 g 2 h 2 , a 3 = %bib 4 g 3 hi, a 4 = — %b\b 3 g 4 h\, as = 
-!*>i b\9 b h\, h 0 = b 4 gr, h x = /t 2 [(l - 2/igro 2 )] -1 and h 2 = [2 - hlw 2 ]~ l . Similarly, 

a R = — 2C m tR* - 4a 2 r 2 (R*ft - flR*) - 8a 3 r 3 (ft 2 R* + R*ft 2 - ^{R*fl 2 }£) 
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( 41 ) 


-16o 4 r 4 (flR*ft 2 - n 2 R*ft) - 32a 5 T 5 {R*ft 2 }(ft 2 - ^{ft 2 }^) 

and the coefficients have the same form as those from the a s solution with the parameter 6j replaced 
with &2 • 

In two dimensions the problem is less complicated since the number of tensor products necessary 
to express the solution is vastly reduced. In this case S,ft,R are two dimensional tensors as 
the mean quantities are 2-D. First it is necessary to recast the tensorial equation in terms of the 
proper traceless 2-D tensors, that is 5* ; ({ it)) = Si-j((u)) — ^5pp((w))^ t -j\ Rij{p-,p) — B, : (p,p) ~ 
±R pp (p,p)6\f . Here, the two dimensional Kronecker symbol is 8 W = [tfj 2) ] = 1 for i = j = 1,2 and 
0 otherwise. The pressure strain model becomes: 


n* - P e*j = -CjP ? a tJ A r F ar + p{k) ^ + ^xM, 2 


£,■(<«»+ 


[1 - c 3 + xMt 


«*£;(<«» + S* p {( u )) a p] ~ 3. S; q ((u))a p<l 6 tJ 
[1 — C 4 — X^ 2 ][ a «p2pj(( u )) — Ilip(( u )) a pj]T 


^xM?S pp ((u))a i: - 


5 3q S pp (( u )) 


3 


^ (2) 
t'n 


2 


2[1 — C3 + X Mt 


x vp 


_ fe 121 ) _ la 16 

3 2 J 3 pq 


II 

3 


8 W 

1 8 

2 


0 


S PP {(u)) 


BrF 0r 


(42) 


It is important to stress the fact the both expressions (2-D and 3-D) of the pressure-strain correlation 
model give the same result when applied to a two dimensional mean flow. The recasting is necessary 
to take advantage of the simplifying properties of this particular case. Also, 


12 12 

M*j = ---T u vo(k)R*j + - -^T u vo(k) R pp 




(2) 

'■J 


(43) 


Using the previous information the final equation becomes: 

-{»£*}«) -MaO-fia) 

'Hi 


a = -gr 


iqS * + b 2 R* + 63 ( aS* + S*a 

’ (8 8 < 2 > 


-(•(s-r 


'6 8 < 2 > 

a 3 _ 2 



(44) 


with h = \-B T F 0r {\ + ^xM?), b 2 = ^1t u j/ 0 ,6 3 = 1 - B T F 0T (l-C 3 + xM?),b 4 = l-B r F r (l- 
C 4 - X M 2 ), 65 = 61 S pp + b 2 R pp , b 6 = b 3 S vp and g having the same expression as in the 3-D algebraic 
equation. The polynomial solution a = £.v C X T X is based on only five independent tensors 


rj\0 


8 8™ 


T 1 = S*, T 2 = S*n - ft{S*}, T 3 = R*, T 4 = R*ft-fiR* (45) 
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and five non-zero independent invariants: 


<r 2 = {S* 2 }, vo 1 = {ft 2 }, {R* 2 }, {R*s*}, {S*K*n}. (46) 

The assertion that there are no other independent tensor products or invariants can be verified 
using the following 2x2 matrix identity: 

2abc = bc{a} -f a{bc} + ac{b} - b{ac} + ab{c} + c{ab) - 

c{a}{b} - a{b}{c} + ({ac}{b} - {acb})£ (2) . (47) 


To obtain the explicit solution to the algebraic equation in the anisotropic Reynolds stress a pro- 
cedure similar to one devised by Pope.'^ We define the matrices 5x5 ^ ^ such that 


T”S* + S*T” - -{T”S*}<5 = ^7f A T A 

3 < 


Jf> 6 W' 

T v I 

l 3 2 


T v n* - n*T v = ^ a t a 

A 

'h <5 (2)A 


2 


-sFU-2 


(48) 


The elements of the matrices are determined from the above equations by making use of matrix rela- 
tions stemming from the Cayley- Hamilton theorem. Next, the coefficients of the tensor polynomial 
are obtained from 


C A = -gr 


61^1 \ 4~ 4~ ~ ^5^oa be 


(49) 


The resulting model for the anisotropy of the Reynolds stress tensor is: 

[ Q 2 S* + (Qi + Qz)bzg f2 T( * 2 — 4" Q-2b 4 g f\T(S H — £1 S ) 

-2C' r [R* + 6 4 tf/,r(R*Q - OR*)] 


a = —2 C^t 


(50) 


The parameters C M and C ' are given by: 



r hgh/2 

l-l(b 3 gra) 2 hf 2 -2(b 4 gf,TZvy 

(51) 


r , b 2 gf i/2 

11 ~ 1 - 2{b 4 gfiT&) 2 

(52) 

and 

b 2 r{S*R*) rt , , {S*R*ft}l 
Qi - 1 + -r — T + 2gUrb 4 X ~ ~ 2 

01 L & ° J 

(53) 
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(54) 


Q 2 = 1 + 


2 (hgraffih (Q _ jx 

31 -2(b 4 gf l Txv) 2 ^ 


36 x 


9f2T 


65 1 - 2(6 4 g/lTCT) 2 
~ 3 _ Ma ZgfiTP 2 


with /1 = (1 + begr/6) 1 and / 2 = (1 - b^gr/d) 


(55) 


2.3.2 ATFM 


The temperature flux transport equation reads: 

dxj 


dp(u" { T") &p(u,T"){u : ) _ d(p(u]u t T ) + pT fr, - T Qjiju )) 

- + ~ ' 


dt 


+p'— - p l + (u J r > ax, J 1 dx, + ax, 




ax, 


a 

ax. 


a „ ar 

• n. 


PriZe *ax, 


_ _P_d±dJT) _ _P_du L dT^ _ a , ( V)^I_ 
PrRe ax, dx., PrRe 0x 3 dx : 3 ' dx } 


-(7 - up + <d («:g ) + (^g )) + 7(7 - ^ 


, a«” 

ax, 


a«*: 


(56) 


Modeling of T" (following the methodology used by Ristorcelli 47 in deriving a model for u” ). From 
the instantaneous the mean is subtracted. 


j t { P T - p(T)) + j^{pT «, - p(T)(u,)) = RHS = -^ (p^ ax, (T T) , 


, / an, du, 
-7(7 -UM 2 


+ 7(7 


1)M 2 + 7Ce(^p - w P ). (57) 


Obviously RH S = 0. After expanding the differences on the left hand side, the terms are expressed 
in non-conservative form and the resulting equation is divided by the instantaneous density. After 
using the approximation 1/p = (1 - p'/p + • • •)/? the equation is averaged and th e terms oforder 
> \f(^/p are eliminated. For example, the difference term {u t T")-u,T'' = p'u.T’/p ~ 0(\fp*/p) 
is discarded. This results in: 


or = T - d A-7 a J21 + 2222. (58) 

Dt ax, J ax, p ax, 


Precedent contributors 49 47 have represented the term correlating fluctuating divergence with veloc- 
ity fluctuation by a linear relaxation model. Because both the divergence-temperature correlation 
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and T” are scalars the following relation is valid on dimensional grounds: 


( 59 ) 


where = M t r is the acoustic time scale. Next, assuming that the quantity T" / yJ{T” 2 ) reaches a 
near-asymptotic state we have 



DT 

Dt {T” 2 ) 


(A-i) 

\P J 


(60) 


where P$ is the production of temperature variance. Combining the the equations 58, 59, 60 the 
following expression is obtained: 


. n &p _—d(T)' 


(61) 


where r t = M t r/[ 1 + Mt(P$/(p e$) - 1)]. Similar to the treatment of the Reynolds stress, 


(T^Tc 

the temperature flux equation is transformed into the equation for the correlation coefficient 

«T") 


t h = 




(62) 


and the transport equation reads 

_D$i 1 

P 


'&T* * 


(k) dTj Ur 2 ) dT ) 

Dt y^r^) \o*i ^ 2 V (T” 2 ) 2 V ( k > d r. 

M + V + p'd — pt, 


( P-d _ d, 


(T 


” 2 \ 


P frf 


2r V P ts 
+ 


1 + 
1 


\fwr*) 


P £s 


(63) 


where the notation D/Dt indicates the convective transport, Tf r Tj and Tj denote turbulent 
transports of the temperature flux, the temperature variance and the kinetic energy, respectively. 
Moreover P# = -yJ(k)(T" 2 )Pjd(T) /dx j is the production of temperature variance and the remain- 
ing quantities are the normalized production, pressure-gradient correlation and flux dissipation: 


P,* = -yJ(k)(T 2 ) 


2 d(T) 1 

<*>K + + ^( 5 0 «“» + « 0 -((“)))+ jO, ■*%«»>) 


(64) 


The pressure-temperature gradient model 

~P fii? = A 4- 2 p!ijk(Sjk((u)) + ft,/t((u)) 


(65) 
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The linear nature of the temperature equations imposes that the integrals Xijk as well as the vector 
be linear in the turbulent flux. Therefore 


Xiik 

\A*><r») 


A x = —C\$p cd, 

= PlSijflk + /?2 (Sik'&j + biktij) + + f 2{ (Lik'd J + dik'd j) + fs^ijdkpdp 


(66) 


The symmetry constraints are satisfied for the above form for To determine the coefficients 

more constraints are applied: the normalization condition which translates into 2„j = ( u"T "). 
The incompressibility condition is replaced with l t kk/\j{k)(T" 2 ) = fM 2 / 2(d, + a^dk), where the 
constant £ has to be calibrated from DNS. These conditions are still insufficient for the deter- 
mination of the constants /i, / 2 , / 3 . These coefficients were chosen so that Iijk proposed by 
SLC89 is recovered in the incompressible limit. Realizability criteria is based now on the tensor 
djk = ({djT ){u k T ) - (u J u fc ){T 2 ))/({u p T )(u p T ) - 2{k)(T 2 )) and the linear function of the in- 
variants of djk Fd = 9/2 — 27djj/2 + 9 djj, I Id is the second invariant of the tensor djk, d 2 - = djidij 
and djj — . 


r(u e T")A e - (T 2 )(A ee -^f) + 2 p Z^u’ 2 ) = CFf, as F D — * 0 


3 

V’2 


^({T 2 )l peqe -{uJ)l epq )-+ 0 a sF D 


(67) 


• Q £ 

- yf(k){T” 2 } ~ ~ Cw dJpjd ld A T F^ r +p |(ci + c 2 )S* J ({u))d^ + (cj - c 2 )fl,j((u))i?j tf 
+(c 3 + c 4 )a,jSjk( (u))dkd + c$ajkS*j{(u))dkti + (c 3 - c 4 )a,jfi jfc ((u))^ 


■b c 5 a jk^»j(( w ))^fct? F c 4d jkSjk({u))di-g + ( t? tt ^ -(- a.ik'flki!) )S pp ( ( u ) ) 


BrFfr 


( 68 ) 


where C lt > = 3.2, c 4 = 4/5 - £M 2 /5, c 2 = -1/5 + 3£M t 2 /10, c 3 = 1/10 + £M?/2, c 4 = -3/10 + 
3^Af 2 /2, c 5 = 1/5- £Af ( 2 . The parameters are a r = 0.1, (3 r = 0.5, A r = min(F5" r , 0.1 -Or ) 
and B r = min(/ 1 £) /?r , 0.1~^ r ). The coupling between the species and temperature reflected in the 
temperature flux equation for the reacting case with heat release is neglected for simplification 
purposes. Using all the present information the equation for d is: 


+ D d Atf + Ctf = 0 


where the coefficient D$ = 2 rh$, with 


h# = 2 A t F%C^ - 1 + (1 + 2c 4 B T Ffr)X = - + - l) + 

(s p c, \p €# J 


Y (1 - xM 2 b t f%)S pp ({u}) + M + p ' d p lc 


P ( s 


-1 


(69) 


(70) 
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where r# = 2{k)(c$/(( s (T" 2 )) is the time scales ratio. The vector term reads: 

= Dtfi 


(k) 


{ T ” 2 > L 


" 1 " 2 ^^* ^ 2 Tt'R-ik 


d(T) 


dx k 


( 71 ) 


Due to the nonsymmetric properties of the second order tensor A an anisotropic turbulent diffusivity 
tensor will be obtained in the final model. 

Aik = X + § TZ, k = [1 - (c, + <*)*„*£]$£((«)) + [1 - (c, - c 2 )F r F^]D, fc ((u)) - 


|(C3 + c 4 )ajj5J*((u)) + c 5 a k ,Sp{(u)) + (c 3 - c 4 )a IJ fl jfc ((u)) - c 5 a A: jn j ,((u)) 

1 

3< 


+^M ( 2 a lfc 5 p * p ((«)) 


+ ^n ik . 


Now, the solution of the system is conveniently represented in the matrix form: 

— -M _1 


(72) 


(73) 


where M denotes the matrix [<5 + D# A]. 


To provide a computationally efficient algorithm, the matrix M is inverted analytically. This is 
achieved via the use of the Cayley-Hamilton theorem and yields an expansion on the minimal 
vectorial basis for this problem: 

2 


W = a nA n C# 

n—0 


(74) 


2.3.3 ASFM 

The methodology used for the temperature flux is now applied to the scalar flux which is transported 
according to the equation: 


dp{u { Y' a ) + dp(u"X)(u 3 ) _ d(p{uyX) + PX^j ~ y>ji («")) 


dt 


dx, 


dxj 


n (uu) d{Ya) I (ur) d{Ui) \ ywda ]t {(u)) 

+p dx, p x u ' ] dx, + {u ^ o) ~d7~ r^dT^^x 


+ 


d 

dx j 


n „ dY, 

-u. 


Sc Re ' dx. 


P du” d(Y a ) p du” dY” „ d\ 

cr,i{u ) 


ScRedx, dx, ScRedxj dx, 


dx, 


where the reaction source term can be approximated by: 

1 1 


u c 7 a = -Da exp 


-Ze 


(T) Tj 


p 2 ({uX)(Yp) + (uX)X) + XXYp)). (76) 


(75) 
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The model derivations leading to expressions for the averaged Favre fluctuations, pressure scalar 
gradient correlation are identical with the ones presented for the temperature related quantities. 
For Y” we have: 



,-V-l 


(77) 


where r a = M t r/[ 1 + e a ) - 1)]. The correlation coefficient 


evolves according to 


Via = 


(«X) 

y/(k)(YZ 2 ) 


(78) 


-DVia 1 (OT* ^ r^T 3T J Via RyP)dTj\ 

Di 2 V OC 2 ) dx 3 2 V <*> dxj 

x + M + V + - p Cc \ + 

pt* /J 

/ [j°io + a ~ P ft a + (79) 

y/(k)(YZ>) 

where the notation is similar to the temperature flux case, for example 

Pa = -y/(k){Y;*)v ja d(Y a )/dz j is the production of scalar variance. The terms that appear extra 
wuth respect with the previous case are due to the reaction source term and they are: S Q — Y ” ) 

which is the chemical source term in the (Y^ 2 ) equation and similarly for the flux equation, 

S; a = -Da - — L =(y,- 0 (Vg) + finO'a) + l.cftJxp))- (80) 

^wk 2 > 


-p 




(O \p fa 


Pa 


- 1 + 


P (a 


+ 


Via 

~2r 



Here = (u" t Y” Yp)/ (k) (Y” 2 )(Yp 2 ) is the normalized covariance flux vector. It is worth 
mentioning that for the pressure scalar gradient model the realizability is based on the tensor 


d ik = {(u]Y' a ){u k Y:) - (u”uI)(Y” 2 ))/((u;Y")(u p Y") - 2<*)<0). The final form of the model is 


the same as for the pressure temperature gradient model, replacing the temperature flux with the 
scalar flux along with the appropriate normalization. This procedure leads to an algebraic system 
of equations for the two unknown vectors <p lQ and 9 ,^. For the mixing case the two vectors are 
uncoupled and the solution is matricially the same as the temperature flux. The reacting can be 
solved exact, but due to chemical coupling the solution will be too complicated to be of practical 
use. Therefore either the mixing solution can be utilized or a perturbation solution of the react- 
ing case, the small parameter t - Ttjji 1 . The fluxes are decomposed as <p xa = + £<p) a and 
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(pip — tpt-p + s<p}p. The solution is obtained from the system of equations 


<p a + D a (A f + el Z)(p a + B a <pp + C Q — 0 
<Pp + Dp( A f + elZ)<pp + Bp^p a + Cp = 0 


( 81 ) 


where the coefficients are 

D a = 


B a — Dap ex p 

And the vector terms read: 

C tQ — D a 


2 rh a 

1 + 2 Dap exp -Ze Th a {Y 0 ) 

1 1 


-Ze 


(T) 


(Y a )D a , 


(82) 


{k) f 2 

(y ”2^ 4" 2^^' — ^2 TtTlik 


d(Y«) 


+ 


Da p exp 


— Ze 


1 

W) 


t d 


dx k 


(83) 


The terms indexed with (3 are obtained from the a indexed terms by the permutations a — ► (3 and 
/? — ► a where necessary. The parameters h a and hp have the same form as h$ (Eq. (70)) with the 
respective change of index. The solution of the system (81) given in matrix form: 


¥>° = -M- 1 [(* + D 0 A)C a - B a C p \ 
V>$ = -M- 1 [(6 + D a A)C 0 - B 0 C O } 


= -M- 1 
= -M - 1 


(6 + DpA)D a K<p° - B a D p n<pi 
(6 + D Q A)D 0 TZ<p$ - B 0 D a n*p% 


(84) 


(85) 


where M denotes the matrix [(1 — B a Bp)6 + ( D a + Dp) A + D Q DpA 2 ]. Using the methods men- 
tioned in the temperature flux section and described below results: 


<P Q = £ «nA n C a + “nA n C 0 . 


(86) 


n=0 


n=0 


2.3.4 Explicit Solution 

The procedure leading to explicit solutions for the turbulent fluxes vector is described. Consider an 
arbitrary three-dimensional second-order tensor Q = [Qtj] and the corresponding Kronecker tensor 
6 = [^j]. According to the Cayley-Hamilton theorem, this matrix satisfies its own characteristic 
polynomial: 

Q 3 - 7qQ 2 + IIqQ - IIIq6 = 0 (87) 
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where / Q = {Q} = Qu, IIq = i[{Q} 2 -{Q 2 }] = \[QnQii -QijQji], IIIq = ±[{Q} 3 -3{Q}{Q 2 } + 
2{Q 3 }] = WQaQiiQ hk - ZQaQjkQkj + 2 QijQjkQki] are the three tensorial invariants. Multiplying 
the characteristic polynomial with Q -1 and solving for the inverse we obtain: 

q~ 1 = 77 hq 2 - 7 <?q +/ w ( g8 ) 


This relation can be used now to find explicit solutions to the problem considered here. We can 
write: 



9o = -(^ + G)- 1 C. 

(89) 

Hence: 

= 777 [G + (2 - I(,+g) g + (1 - h+G + Hs+g)^]- 

Mh+G 


05 + G)- 1 : 

(90) 

It is easy to show that I$+g 

= lG + {b}- Ih+G ~ 2 Iq + IIq + {S}, Uh+G = Ig+Hg + IHg + ^- 

Therefore the normalized turbulent flux vector takes the form: 



*ra — «oC + ajGC + a2G 2 C 

(91) 

with the coefficients: 

1 + ^G + IIg 

a °~ 1 + I g + IIgF I IIg 

(92) 


1 + /g 

ai “ 1 + / G + // G + 1 IIg 

(93) 


1 

° 2 — — l + Ig + Hg + m G 

(94) 


The reacting case is somewhat more complex. Nevertheless, by following the same procedure 
explicit solutions are obtained: 


ifc* = GoC 0 + Uq Op + ci\ AC 0 + a\ACp + o,2A^C Q + a 2 A 2 C/3 ( 9 b) 

<P0 = b 0 C a + b' o C 0 + bi AC a + b\AC p + b 2 A 2 C q + b' 2 A 2 C 0 , (90) 

with the coefficients: 


Go — ~ 


F 0 (F a + + B a B 0 [yft(D 0 F o - ED a )D 0 - E(E + i^±D a D 0 ) 

(1 -B Q B 0 )(F a F 0 -E 2 B a B 0 ) 


(97) 


Gq — Bq f 


F a {F 0 + D 2 0 ^p - ) + D a yp(D o F 0 - ED 0 ) - EB a B 0 (E + ^D O D 0 ) 
(1 - B a B 0 )(F a F 0 - E 2 B a B 0 ) 


(98) 
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( 99 ) 


Bg Bjj [E[Dp - Dp) j- FgDp\ - DgFp 
D a (F a Fp - E 2 B a Bp) 


fl i - d 

^ a 


B a DgFj3 — DpFg - E(Dj) — DJUb) 
F a Fp - E 2 B a Bp 


( 100 ) 


«2 = 


D a Fp - EDpB a Bp 


with the shorthand notations: 


E a = (1 - B a Bp)(D a 


Fp = (1 - B 0 Bp){Dp 


E a Fp - 

E 2 B a Bp 

D D a Fp - EDp 

" a F a F 0 

- E 2 B a Bp’’ 

, {A 2 } 

1 B a Bp 

° 2 

D a Dp 

, {A 2 } 

1/5 — 

1 B a Bp 


)~Bl 


{A 3 } 


-D/3 F) a 


)2 {A 3 } 

^ 3 


( 101 ) 


( 102 ) 


(103) 


(104) 


1 1 /A 3 ! 

^ = (^ + ^)(1 - B a B p ) - D Q DpE-l. 


The coefficients £>, are obtained from the a,’s through the permutations a 
a' 0 — *• b 0 , ai — * b\, a\ -*■ b u a 2 -* b' 2 and a' 2 — ► b 2 - 


(105) 

(3, /3 -> a, a 0 -+ b' 0 , 


2.4 Results and Conclusion 

The theory built in previous chapters will be used to simulate numerically a non-premixed, tur- 
bulent reacting with heat release, spatially developing mixing layer over a wide range of Mach 
numbers. The simulations are performed on a uniform grid in the computational space. By means 
of a coordinate transformation the mesh is transversally compressed in the physical space in the 
region corresponding to the actual mixing layer and the equations are solved in vector form. The 
numerical solution procedure for the integration of the governing equations is based on a Gottlieb- 
Turkel predictor-corrector finite difference scheme. 50 The method has dissipative properties and it 
is second-order accurate in time and fourth-order accurate in space. The interest at the present 
time lies in the steady-state solution and therefore, to accelerate the convergence towards steadi- 
ness, a local time stepping technique is used. The convergence criteria is imposed to be that the 
reduction of the steady-state residua] in average sense attains a minimum acceptable level. More 
specifically, the simulation is considered at steady-state when is observed at. least a 1.5 order of 
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magnitude reduction of a global quantity, such as the absolute value of the residual averaged over 
the whole domain. Figure 2.1 shows a typical evolution for the the previously mentioned quantity 
obtained with the Gottlieb-Turkel scheme. Although the criteria is quite stringent it is known 51 
that predictor-corrector type of schemes are not able to achieve better rates of residual reduction. 
The set of initial conditions is obtained by propagating the inflow conditions throughout the entire 
domain, hence, in this procedure, the flow has to sweep at least one time the domain to obtain a 
meaningful result. The boundary conditions are set according to the elliptic nature of the problem 
on all four boundaries. The inflow BC specifies smoothed step or smoothed hat profiles for the 
primary variables. At the outflow and outer boundaries zero gradient boundary conditions are 
applied for their nonreflective properties in relation with the outgoing waves. The grid overlay- 
ing the computational domain of 120^ x 60(5^ had 128 x 64 points, where the vorticity thickness 
= (ui - U 2 )/(du/dy) max . In the first stage of this work we have concentrated mainly on hy- 
drodynamics. Test simulations were performed with simplified versions of the algebraic model for 
the Reynolds stresses. The particular set of conditions for the two streams of air were a prescribed 
velocity ratio r v = 1/4 and equal thermodynamic properties. The issue of reduced spreading rate 
of the shear layers with increasing free stream Mach number is well established experimentally. 
The compressibility effects parameter called convective Mach number correlates with the growth 
rate normalized by its incompressible value at the same ratios for velocity and density 52 and had 
the expression in this case M c — M (1 - r v )/ 2. The fully developed shear layer at high Reynolds 
numbers grows linearly and the spreading rate can be expressed as d6 u /dx = C$( 1 — r v )/(\ + r v ) 
where 6 = 6(x) is the thickness of the shear layer based on the normalized velocity profile and C$ 
is constant (approximatively). Figure 2.2 represents the downstream-evolution of the shear-layer 
width. The linear growth is attained after a phase of development near the inlet. The confirmation 
that a fully developed regime with linear growth has been installed in the flow, stems from the 
self-similar property of the velocity profiles and normalized Reynolds stresses. As it can be seen 
from figures 2.3, 2.4, 2.5 the profiles, plotted in similarity coordinates, collapse for each axial co- 
ordinate considered in the outflow region of the domain. The temperature profile shows the same 
self-similar behavior as the other mentioned quantities and as a result of the velocity gradients 
displays an increase in the middle of the layer. The correlation profile of the normalized C$ versus 
M c has yet to be completed. The simulations will be continued with the full set of algebraic models 
and tests will be made for the entire range of interest of convective Mach number. 

The purpose of the efforts in this part of our activities was to develop closures for the “second order 
moments” in the contexts of both RANS and LES of high speed turbulent flows. In particular, 
we have obtained a complete set of explicit algebraic models derived from a hierarchy of second- 
order moment closures that are valid for compressible turbulent flows. The primary methodology 
based on the Cayley- Hamilton theorem and its corrolaries was developed during the past twenty 
years 18,19,1 and the present work extends it further. The models are based on new compressible 
closures for the pressure-strain and pressure-scalar gradient correlations. Explicit algebraic relations 
are provided for the Reynolds stresses, velocity-temperature and velocity-scalar correlations in both 
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non-reacting and reacting flows with heat release. As theoretical novelty we mention that density 
gradient, pressure gradient and mean dilatation effects are included in the models. Also, the role 
of the turbulent Mach number and Damkohlernumber is exhibited as well as compressibility and 
variable inertia effects for application of the models to turbulent flows with nonpremixed reactants. 


3 Monte Carlo Large Eddy Simulation of Reacting Turbulent 
Flows 

3.1 Introduction 

The purpose of the efforts in this part of our activities is to develop and implement a robust 
computational procedure for LES of turbulent reactive flows. The procedure is based on a Monte 
Carlo solver for the Probability Density Function (PDF) of Subgrid Scale (SGS) of reactive species. 

The purpose of the efforts described in this report is to develop and implement a robust compu- 
tational procedure for Large Eddy Simulations (LES) capable of capturing the intricate physics 
associated with turbulent reactive flowfields. LES is considered somewhere between Direct Numer- 
ical Simulation (DNS) and Reynolds Averaged Navier-Stokes (RANS) computation . 53-58,11,59 Over 
the past thirty years since the early work of Smagorinsky 60 there has been relatively little effort, 
compared to that in RANS calculations, to make full use of LES for engineering applications. The 
most prominent model has been the Smagorinsky eddy viscosity based closure which relates the 
unknown subgrid scale (SGS) Reynolds stresses to the local large scale rate of flow strain. This 
viscosity is aimed to provide the role of mimicking the dissipative behavior of the unresolved small 
scales. The extensions to ‘dynamic’ models 61-64 has shown some improvements. This is particu- 
larly the case in transitional flow simulations where the dynamic evolutions of the empirical model 
‘constant’ result in (somewhat) better predictions of the large scale flow features. 

A survey of combustion literature reveals relatively little work in LES of chemically reacting turbu- 
lent flows . 55,65 It appears that Schumann 66 was one of the first to conduct LES of a reacting flow. 
However, the assumption made in his work simply to “neglect” the contribution of the SGS scalar 
fluctuations to the filtered reaction rate is debatable. The importance of such fluctuations is well 
recognized in RANS of reacting flows in both combustion 12 ’ 15,11,17 and chemical engineering 67-70 
problems. Therefore, it is natural to expect that these fluctuations will also have a. significant 
influence in LES. 

Modeling of scalar fluctuations in RANS has been the study of intense investigations since the 
pioneering work of Toor . 71 The aim of statistical moment methods is to provide a closure for 
these correlations in terms of mean flow variables. Because of the lack of models wdth universal 
applicability to accurately predict the scalar correlations in turbulent reactive flows, simulations 
involving turbulent combustion are often met with a degree of skepticism. Another approach which 
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has proven particularly useful is based on the Probability Density Function (PDF) or joint PDF 
of scalar quantities . 72, 55,73 This approach offers the advantage that all the statistical information 
pertaining to the scalar field is embedded within the PDF. Because of their capabilities, PDF 
methods have been widely used in RANS for a variety of reacting systems (see Dopazo 73 for a recent 
review). A systematic approach for determining the PDF is by means of solving the transport 
equation governing its evolution . 74 In this equation the effects of chemical reaction appear in a 
closed form. However, modeling is needed to account for transport of the PDF in the domain of the 
random variables. In addition, there is an extra dimensionality associated with the composition 
domain which must be treated. An alternative approach is based on assumed methods. In these 
methods the shape of the PDF is assumed a priori usually in terms of low order moments of the 
random variable(s). Obviously, this method is ad hoc but it offers more flexibility than the first 
approach. Presently the use of assumed methods in R ANS is justified in cases where there is strong 
evidence that the PDF assumes a particular distribution . 75-77 

Despite the demonstrated capabilities of PDF methods in RANS, their use in LES is limited . 55,78,59 
The first application of PDF-LES is due to Madnia and Givi 78 in which the assumed Pearson family 
of PDF’s are used for modeling of the SGS reactant conversion rate in homogenous flows under 
chemical equilibrium conditions. This very same procedure was also used by CR 79 in the LES of a 
similar flow. The extension of the model for LES of nonequilibrium reacting flows is reported by 
Frankel et a /. 3 for LES of reacting shear flows. While the generated results are encouraging, they 
do point out the drawbacks of assumed PDF methods. These can be overcome only by considering 
the PDF transport directly. 

The approach advocated here is to solve the transport equation for the PDF. Because of the added 
dimensionality due to the compositional variables, solution of the PDF transport equation by 
conventional numerical methods is possible in only the simplest of cases . 80 An analysis performed 
by Pope 81 suggests that the solution of the joint velocity-scalar PDF by finite difference methods 
is impractical for more that three scalars. 

The numerical solution of the subgrid PDF may be accomplished by means of a u Monte Carlo” 
scheme. The use of such schemes in RANS has proven very effective , 72 however no attempt has 
ever been made to utilize Monte Carlo schemes in the context of LES. Two classes of Monte Carlo 
schemes exist. In the Eulerian type scheme, the PDF within the subgrid is represented by an 
ensemble of M computational elements at fixed grid points. These elements are “transported” in 
physical space by the combined actions of large scale convection and diffusion (molecular and sub- 
grid eddy). In addition, transport in compositional space occurs due to the processes of chemical 
reaction and molecular mixing. Preliminary implementation of an Eulerian Monte Carlo method 
for LES of a non-reacting mixing layer has been performed in unpublished work. The Smagorinsky 
closure was used for the modeling of the subgrid eddy diffusion and the stochastic model of Curl 82 
was utilized for modeling of the molecular mixing. Unfortunately, the results were quite discour- 
aging. The major difficulty with this formulation lies in the numerical implementation of the large 


25 



scale convection. Due to the nature of the grid based scheme, excessive artificial diffusion is created 
which greatly degrades the solution of the large scale structures. It is important to realize that 
the errors induced by this scheme are not at all due to the PDF formulation itself but rather to 
the numerical implementation of the closed mean convection term. A remedy for the problem is to 
divorce from the Eulerian discretization and to invoke the Monte Carlo solver for the LES-PDF in 
a “grid free” Lagrangian manner. 

In this work we provide a computational methodology for solution of the PDF of SGS scalar 
variables in LES of reacting flows under nonequilibrium chemical conditions. The solution procedure 
involves the transport of N Lagrangian elements within the “whole” computational domain of 
interest. The advantages of Lagrangian numerical procedures in reducing numerical diffusion in 
DNS are well-recognized. 83-86 In this Lagrangian framework, the elements are free to move anywhere 
within the domain. The particles carry information pertaining to the scalar field only; the LES 
of hydrodynamic variables is conducted by conventional Eulerian finite difference procedures. The 
effects of convection and diffusion are to move the elements in physical space, while the effects of 
mixing and reaction are to modify the compositional makeup of the elements. 


3.2 Governing Equations 


The primary independent transport variables in a compressible, two-dimensional flow undergoing 
chemical reaction are the density p, the velocity vector u t , the total specific energy £, the pressure p, 
the temperature T, and the species mass fractions f Q (a = 1,2, ..., N s ). The conservation equations 
governing these variables are the continuity, momentum and species mass fraction equations, along 
with an equation of state relating thermodynamic variables. They are expressed as: 


Continuity : 

dp djm 
dt dx, 

Conservation of momentum : 

dpuj dpu t Uj _ dr tJ 

dt dx, dx. 

Conservation of total energy : 

dpE dpu x E drijUj dq t 

dt + dx, dxi dx-i 

Conservation of chemical species : 

dpfg dpuja _ch/£ , - 

dt dx, dx, ° 


(106) 

(107) 


(108) 


(109) 
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Equation of state : 


N s 

p = pR°Tj2f>/M, 

1 = 1 


(110) 


The total specific energy is given by : 


A ’ a n v 2 + V 2 

1=1 ” 


( 111 ) 


and the enthalpy of species i is defined as 

h t = /i° + [ T c p ,{T')dT' 
JT r 


(112) 


Additionally, the viscous stress tensor r,y, heat flux q{ and mass flux «7 t ^ of chemical species a are 
given as: 

rhi rhi r) / n i. 

(113) 


, , dui duj duk 

r u -^p + p ( — + ^- ) + ^A_ 


9i ^ 


&r, 




(114) 

(115) 


3.2.1 Modeling of Unresolved Scales 

The aero-thermodynamic equations of the previous section constitute a complete set of governing 
equations. Unfortunately, due to the limited power of today’s computers, it is impossible to accu- 
rately solve these equations for typical engineering problems. The great variation in length scales 
would require grid resolutions that would be too prohibitive for even the fastest of today’s super- 
computers. RANS provides the engineer with an alternative; instead of obtaining a fully resolved 
solution which can be afforded in only limited cases, time averaged solutions which do not attempt 
to resolve the the fine structure of the turbulence could be attempted. Due to the non-linear nature 
of the equations, the time averaging procedure yields unclosed terms which have been the focus of 
much attention in the past. On the other hand it may be desirable to resolve some of the lower 
frequency turbulent structures. Instead of averaging over all time (and implicitly length) scales, 
LES attempts to resolve the larger, energy containing eddies. Because only the finer turbulent 
structure is modeled, it is expected that LES models would be more universal in application. This 
is accomplished by use of a local spatial filter: 

/ + OG 

0(x, t)G(x. f - x)dx' (116) 

-oc 
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The fluctuation about the filtered value is given by 


<t>' = 4> ~ 4>- 


(117) 


The filter G(x) can take many forms. In this work we have elected to work with a local volume 
average box filter. For compressible flows, it is desirable to work with Favre averages: 


4>=^- 

p 


(118) 


and the fluctuation about this mean is denoted by 

<p" — <p — <f>. ( 119 ) 

For compressible flows with reaction it is convenient to work with the density weighted mass fraction 

Fv = Pfa/P ( 12 °) 


so that T 0 - f a . Note that while 


N. 

fa — 1 5 

k = 1 


A'.- 


Y^ :Fa = P/P / 1- 


(121) 


( 122 ) 




When the LES Favre averaging procedure is applied to the governing transport equations, the result 

(123) 


is: 


dp (Tim _ 

dt dx t 

djntj dpu t uj _ dj\j _ dT,j 
di + dxi ~ dx, dx i 

dpE dpu,E _ dT tl itj _ d % _ dQi 

dt. + dxi ~ dx i dx, dx i 

dp Jo dpuja _ _dJ]_ _ dM? £. 
dt dxi dx, dx, 


(124) 

(125) 

( 126 ) 


where 


T,j = p(u,Uj - u,Uj) (127) 

Q, - p(u,E - UiE) (128) 

M° = p(ui7a ~ Vila) ( 12 °) 

are unclosed terms and therefore a model must be provided to account for their effects. 

In this preliminary work, which is restricted to no heat release and low compressibility, the in- 
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compressible Smagorinsky eddy viscosity subgrid model is employed. The use of an incompressible 
model is justified in lieu of the fact that the density variations are expected to be quite small. Using 
this formulation, the SGS stresses are given by: 

Tij - (6ij/Z)T k k = -2/itSij. ( 13 °) 

where S Jt is the large scale strain rate tensor. A similar eddy viscosity formulation is used to close 
the SGS heat and mass fluxes: 


fi t c v dT 

(131) 

Pr t dxi ’ 


Pi dfa 

(132) 

Sc t dx t 



The Smagorinsky eddy viscosity is given by: 

Pi = pC 5 A 2 |5|. O 33 ) 

where \S\ = and A is the filter width. In this work, the constants C s , Pr t and Sc t are 

set to the values 0.010,0.7 and 0.7, respectively. 

Thus far we have not yet addressed the issue of how to deal with SGS scalar correlations in the 
filtered chemical source terms. While the SGS terms discussed in this section are of a convective 
nature and they can be reasonably well modeled by a diffusive process, the same cannot be said 
for the unclosed terms in the species production rates. Because the physical mechanism of the 
SGS stresses and fluxes is inherently different from the scalar correlations in the source terms, it 
is expected that the models will differ. In fact, when eddy viscosity concepts are extended to treat 
chemical source terms, the resulting models (“eddy break up models”) perform mediocre at best. 
The focus of the following sections is to discuss how the methodology of LES via PDF (hereinafter 
refereed to as LES-PDF) is used to overcome the closure problem of the chemical source terms and 
to develop robust numerical methods for the simulation of turbulent reactive flows. 

3.3 PDF Methodology 

The most common approach to turbulent reactive flow problems in the past has been to solve 
the governing transport equations for the Favre averaged flow variables. As a consequence of the 
averaging process, unclosed terms appear and must be modeled. This type of methodology is 
referred to as moment closure methods as closures must be provided for the unknown correlations. 
An alternate approach is consider the joint PDF of scalar quantities rather than to directly solve 
the scalar transport equations directly. Once the joint PDF is known, all statistical quantities 
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involving the scalar variables can be determined. The average species production rate u 0 can then 
be determined since it is dependent on such scalar correlations. The natural starting point for the 
consideration of PDF methods is to examine the transport equation governing its evolution. First, 
however, it is convenient to discuss some preliminaries. The following is intended only as a brief 
review and is by no means complete. For further information the reader is encouraged to consult 
one of the many fine texts on stochastic analysis such as those by Schuss 87 and Papoulis 88 and 
BiUingsly 89 An excellent reference of PDF methods in the application to turbulent reactive flows 

is provided by Pope. 72 

3.3.1 The Probability Distribution Function 

Let <p{x,t) be a random variable. The possible values that can be assumed by <p constitute the 
sample space. In general the sample space may consist of all of the real numbers, but further 
physical restrictions may restrict allowable values to a subset of the real numbers. For example, 
the thermodynamic temperature T can only take on non-negative values. If we regard 7 as a 
random variable, the sample space of T consists of all the non-negative real numbers. In a similar 
manner the sample space for the mass fraction of any species consists of all the real numbers from 

0 to 1. 

The probability distribution function (also commonly referred to as the Cumulative Distribution 
Function or CDF) is defined for a continuous random variable <f> by: 

W) = P{<P < 0}> (134) 

where P{A} is the probability that the event A occurs. The probability of an impossible event is 
zero, while the probability of a certain event is unity. Some fundamental properties of CDF’s are 


i>(- oo) = p{4> < -oo) = 0, 

(135) 

J>(+ 00) = P{4>< +00} = 1, 

(136) 

B f. * > 0. 

(137) 


The last requirement merely states that the CDF is a non-decreasing function. 

The probability distribution function is useful because we can use it to determine the likelihood 
that a random variable will fall between two values. For example, the probability that the random 
variable <j> will fall between the values 4> a and i/’fc is given by 

P{4' a <4>< Vb) = F+Wa) - -W’fc)- ( 138 ) 
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3.3.2 The Probability Density Function 


The PDF is defined to be the derivative of the CDF: 


U(i') = 


dW) 

di\) 


( 139 ) 


The PDF, like the CDF, is useful for determining the probability that a random variable falls 
between two values. Upon integrating Eq. (139) between ip a and 4’b the following is obtained: 


/ +oo 

. 

-oo 




(140) 


Over a small interval f^Aip is an approximation to P{V’ < <t> < ip + Axp} = A F# and in the 
limit represents the differential probability dF# that the random variable $ lies in the 

infinitesimal interval of width dip in the vicinity of ip. 

Some fundamental properties of the PDF are 


UW) > 0, 

(141) 

U(-oo) = 0, 

(142) 

Ui+oo) = 0, 

(143) 

-f oo 


U(i')di' - i. 

(144) 


— OO 


3.3.3 Determination of Mean Values from the PDF 

Let A{(j>) be a function of the random variable <+>. In general, the expected or mean value of A (<)> ), 
denoted by A(<f>) or </!(<£)> is given by the expression: 

/ +oo 

MWUWdil’. (145) 

‘OO 

^From this equation the PDF can be regarded as a normalized weighing function that assigns 
various weights to each possible outcome in the sample space. The integration serves to sum over 
each of these weighted possibilities to yield the expected value. The term ensemble average is also 
used to denote the mean. 

In particular, the mean value of the random variable <f> is simply 

_ r + OO 

4>= tU{i')di\ (146) 
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and the n th central moment is given by 


/ +oo _ 

W-W)#* (147) 

-oo 

For example, the second central moment, or variance is given by 

/ +°° — „ 

(148) 

-OO 

The variance represents the mean value of the square of the fluctuation and is useful in determining 
how a random variable fluctuates about the mean. A random variable with a large variance deviates 
more about the mean compared to a random variable with a smaller variance. The square root of 
the variance is called the standard deviation and is denoted by the greek letter a. 

In PDF methods, it is often useful to express the PDF a.s the mean value of the delta function: 

U(t)= -<»)>■ (149) 

The delta function 6(ib - <j>) is commonly referred to as the fine grained PDF. Each fine grained 
PDF 8(tp - <t>) corresponds to a realization from the sample space with the fine grain value <f> = ip. 
This utilization of the delta, function is important when we extend the concept of PDF’s to spatially 
filtered quantities. 


3.3.4 Functions of More Than One Random Variable 

Up to the present point functions of only one random variable have been addressed. The corre- 
sponding PDF’s are referred to as univariate since they only depend on one random variable. Many 
situations, however, require the use of functions that depend on more than one random variable. 
The distributions associated with such random variables are referred to as joint or multivariate. For 
example, the reaction rate of a chemical species generally depends upon multiple scalar quantities 
such as several species mass fractions T a (a =1,2, ..., N s ) and the temperature T. Denote <f> n = T n 
for n - 1,2, ..., N 3 and <pm = T where = N s + 1 so that <f> is a random vector of dimension Nj,. 
The reaction rate of species a, is 

d> Q = d> 0 (</>) = W o (<£i,02,...,<AjV*). ( 15 °) 

Notice that the sample space is a dimensional hypervolume. Each point within this hypervolume 
represents a possible composition. Some of these outcomes are more likely than others while others 
may not be possible at all. A multivariate PDF is useful for characterizing the entire statistical 
behavior of the random field. To introduce the concept of a multivariate distribution, we begin 
with a bivariate distribution which involves two random variables. Extension to higher dimension 
PDF’s is inductive. 
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Let <j > i and 0 2 be two independent random variables. The distribution function F # is 
given by 

F+itoWhi’i) = Ptyi < < ^ 2 }, ( 151 ) 

and the bivariate PDF is given by 


/<Ma(lh>lk) = ( 152 ) 

The univariate PDF of <j> 1 may be obtained by integrating the multivariate PDF once: 


/* + 0 © 

f<t>M 1 )= / /*i*a(0i,lfe)<*02. 

J — OO 

(153) 

Similarly, 

r f 00 

f<h (i'l) = / f<h 4>i ( V’1 , v>2 )dip 1 . 

(154) 

Some additional properties of bivariate CDF’s and PDF’s are: 


1 ) = F <t> i^ 2 (V J i,+oo), 

(155) 

-^>2 ( $2 ( V*2 )i 

(156) 

^(- 00 ,^ 2 ) = r^ 10J (0i,-oo) = 0 , 

(157) 

d* 0 ] ( "1” 1 ) 1* 

(158) 

r+oo r -(•OO 

/ / = ii 

J — OQ j — OO 

(159) 

Ui'hWufa) > 0. 

(160) 


Consider a function of two random variables A(4>\,<j>2)- The mean value of this function is obtained 
by integrating over the two dimensional region: 

A = f [ ( 161 ) 

J — OO J — OO 

The covariance C # ^ of two random variables is defined as 

C^ 2 = <(0i (162) 

Higher order joint moments of the product 4>\4> k 2 are given by 

m jk = <4>\4>2>- ( 163 ) 
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and higher order joint central moments are given by 


Vjk = <(4>\ ~ y {4*2 ~ </ > 2) fc >- (164) 

Conditional statistics are important in situations where one random variable takes on a prescribed 
value. For example, the conditional density /<*>, |^ 2 (V’i| 02) given by 

(165) 

is the probability density function of the random variable <f>\ for a given value of the random variable 
<f> 2 . Conditional mean values are given by the expression 

A(4> i|0 2 )= / -4 (^i)/0 1 ^ 2 (V’i|<A2)^V’i (166) 

J — OO 

and is essentially the average of the function A at a fixed value of <j> 2 . 

Extensions to more than 2 random variables is inductive. For example, the mean of the reaction 
rate requires knowledge of the joint PDF 

/ +oc y+oo y+oc 

•••/ / (167) 

-OO J — OO J —OQ 

(When there is no fear of ambiguity as to what random variables are involved, the subscripts for 
PDF’s and CDF’s may be dropped). 

The extension of Eq. (149) to more than one random variable is: 

- 0)> = <<*>( - <7>1 )^( V’2 - 4> 2 ) • • - <t>N<,)>- (168) 


3.3.5 Large Eddy PDF 

The PDF’s considered so far are not fully suitable for LES since they do not contain any information 
about the filter G(x). A “Large Eddy PDF” that is consistent with spatially filtered quantities may 
be defined as follows: 90 

/ + OO 

6{4> - i/>(x, t )]G(x' - x)dx'. (169) 

-OO 

The Large Eddy PDF is therefore seen to be the spatially filtered value of a delta function in contrast 
to the “conventional” PDF which can be defined as the ensemble average of a delta function. In 
either case the delta function is the fine grained PDF introduced in section 3.3.3 

Evaluation of spatial averages and moments is achieved by integrating with the Large Eddy PDF 
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just as ensemble averages and moments are evaluated using “conventional” PDF’s. Consider the 
function A(<f>) of the random vector </>(x,/). The filtered variable A(x y t) is then given by the 
expression: 

/ +oo r -foe 

A(x',t)G(x' - x)dx' = / (170) 

-OO J — oo 


For compressible LES, it is useful to define the Favre PDF fi = pfij < p >. The Favre filtered 
variable A(x,t) is then given by the expression: 


A(x 


/ + CO 

-oo 


A(4>)f L (ip,x,t)d t/>. 


(171) 


3.3.6 PDF Transport Equation 


The transport equation for the joint compositional PDF in a turbulent reactive flow is given by: 


.72 


d[/?(-0)A] d[p(rj>)uJ L ) d[p[ij>)u} q(-0)/l] 
dt 8r t dtp a 


d 

di’o 




djpWJ^l M/l] 

dx{ 


(172) 


Summation over a - is implied. In the context of Large Eddy Simulation, the PDF fa 

has a slightly different meaning than that intended by Pope; the PDF here is the Favre filtered 
Large Eddy PDF as described in section 3.3.5. The first two terms on the left hand side of the 
equation represent convection of the PDF by the mean flow in physical space. The last term 
on the left hand side is due to chemical reaction. Note that this term is closed and requires no 
modeling. This is the major advantage PDF methods have over other approaches. Also note that 
the derivative is in compositional space rather than physical space. This is reflected by the fact 
that the chemical reaction serves to change the compositional makeup of the mixture rather than 
to provide a mechanism for motion in physical space. The first term on the right hand side is due 
to molecular mixing. Molecular mixing, like reaction, provides transport in compositional rather 
than physical space. In general, the mixing term tends to homogenize the fluid and hence lowers 
the scalar variance. The remaining term represents turbulent transport of the PDF by the small 
scales. 


3.3.7 Modeling of Unclosed Terms 

The two terms on the right side of Eq. (172) are unclosed and therefore a model must be provided 
to properly account for the effects they have on the larger scales. In the present work, a simple 
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gradient-diffusion model has been employed for the turbulent transport term: 


p<u"\i]»J L 


.oh 

* dxi 


(173) 


where = pt/ Sc t . 

The modeling of the molecular mixing term has been the focus of intense investigation in the 
past. 82,91,80,92 Many mixing closures fall under the general category of Coalescence/Dispersion 
(C/D) models as characterized by Pope. 93 Curl’s model and the modified Curl model of Janicka 
et al. s0 fall in this category. In the present work, Dopazo’s deterministic relax to mean model has 
been utilized. Kosalv 94 has shown that Dopazo’s model belongs to the general class of C/D closures 
under certain limiting conditions. This model has been selected due to its ease of implementation 
into the numerical solution and its efficiency. Although the use of different mixing models results 
in different behavior of moments of the third order and higher in 0 D turbulent mixing simulations, 
it is expected that the difference under more realistic conditions will have little effect on mean flow 
quantities. 


3.3.8 Langevin Equation 


The basis of the numerical scheme used for the solution of the PDF transport equation relies upon 
the principle of equivalent systems. 72 Two systems that display dissimilar behavior may actually 
have identical statistics. In the Lagrangian solution technique used to solve the compositional PDF 
equation, Monte Carlo particles are distributed throughout the flowfield. Each of these particles 
carries information about the scalar field. Additionally each of these particles obey certain equations 
which govern its motion in three dimensional space. It is important to recognize that the Monte 
Carlo particles are not fluid elements. In fact, while fluid particles follow smooth trajectories, 
the Monte Carlo particles follow trajectories which are continuous but not differentiable. The 
importance, however, of the Monte Carlo particles is that they are developed in such a way such 
that they evolve with the same PDF associated with genuine fluid particles. 

The Monte Carlo particles undergo motion in three dimensional space by convection due to the 
mean flow velocity and diffusion due to molecular and turbulent viscosities. This type of general 
diffusion process is represented by a Langevin equation: 72 


dX,(t) = 


r idr ( i 


[ 2F< 1 

«, + 

P <J%i J 

dt + 

. P - 


1/2 


dW t , 


(174) 


where X t is the Lagrangian position of the Monte Carlo particle. The stochastic term IF, is the 
stochastic Weiner process. The Weiner process is best understood by considering the function 
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Wd{tn) which changes value at discrete time intervals: 72 

N 

WWn) = (At) 1/2 ££i, (175) 

t = l 

where f n (n = l,2,...,A r ) are N independent normalized gaussian random variables and the time 
interval from t = 0 to t = T is divided into N equal subintervals of duration At = T/N . Consider 
the increment 

AHW„_,) = UAi) 1/2 - (176) 

The Weiner process can be defined as Eqs. (175)-(176) in the limit At — ► 0. Note that although 
the process is continuous, it is not differentiable since AWd/At is undefined as At vanishes. 

It must be emphasized that although the Langevin equation given by Eq. (174) is stochastic, Eq. 
(172) which governs the transport of the joint scalar PDF is deterministic. 


3.4 Numerical Solution 

Because of the added dimensionality the compositional variables present in the PDF Transport 
equation its solution by conventional finite difference or finite volume methods is intractable for 
engineering problems. Instead, a Lagrangian Monte Carlo solution algorithm is utilized. It is an 
established fact that while the work required by finite difference schemes increases exponentially 
with added dimensionality, the work associated with Monte Carlo schemes only increases linearly. 93 
Thus Monte Carlo methods provide an attractive technique to solve problems with a large number of 
independent variables. The essentials of Lagrangian Monte Carlo schemes in relevance to turbulent 
flows are due to. 72 

In the solution procedure numerous particles are distributed throughout the domain. Each of these 
particles carries information about the compositional makeup of the fluid. Although Monte Carlo 
particles and fluid particles are fundamentally different their PDF’s are identical. Thus a solution to 
the PDF transport equation can be can be attained indirectly by solving for the spatial location and 
compositional makeup of the Monte Carlo particles. The equation governing the spatial location of 
the particles is the Langevin equation (Eq. (174)) and the processes of mixing and reaction govern 
the compositional evolution. Since the compositional PDF provides no information about the 
density or velocity fields, this information must be determined by alternative means. Conventional 
finite difference schemes are used to solve the governing transport equations for these variables; 
the Monte Carlo procedure is used only to determine the species mass fractions. Additionally, in 
the case of temperature dependent chemical kinetics, the energy (or enthalpy) is solved using the 
Monte Carlo procedure. In the present research temperature independent kinetics are used and the 
total energy equation is solved using finite differences in the usual manner. 
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3.4.1 Finite Difference Procedure 


The methodology used to solve the large scale continuity, momentum and energy equations is the 
fourth order spatially accurate finite difference scheme as developed by Gottlieb and Turkell. 50 The 
generalized transport equation may be written in the vector form: 


r)U OF dG 
dt dx dy 


where U is the vector of dependent conserved variables: 


U = < 


P 

pu 

pv 

pE 


(177) 


(178) 


The vectors F and G are the flux vectors in the x and y directions respectively: 


F = ( 


(pE 


pu 

pu U ^~XX “h 

puv T xy + Exy 

r xy V ~j“ Q x Qz 


(179) 


(180) 


pv 

pv XI Tyx “j“ Eyx 

PVV ~ Tyy + Tyy 

( ( pE — Tyy)l£ — Ty X U (}y — Q y j 

The source vector H = 0 since the species equations are not solved by the finite difference method 
and the body force vector is neglected. 


The Gottlieb-Turkel scheme is a higher order accurate variant of the well known 95 predictor- 
corrector method. For Eq. (177) it is implemented in unsplit form as: 


At 


At 


+ SRW - 7F M - + 8G w+. - roy + mh;, (isi) 

= u '^ ■ “ 8F '-U + 7F U - - 8G?J-| + ? G ',] + A1H-J (182) 


1 


u n+1 = -[U n + U**l 


(183) 


The CFL condition for this scheme requires that the CFL number should be less than 2/3 for 
numerical stability. 
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3.4.2 Monte Carlo Particle Method 


The species mass fractions are determined using the Monte Carlo particle method. A two stage 
Runge-Kutta scheme is used to solve the Langevin equation to determine the particle positions. 
For the general diffusion process governed by: 

dXi(t) = Di(X(t),t)dt + B(X(t),t)dWi(t) (184) 

a Runge-Kutta scheme can be written as: 


X- = Xi + D?At + B n (At) 1/2 C 

(185) 

X" = X’ + D’At + 2T(A/) 1/2 e 

(186) 

A -n+l = I [A? + XT*] 

(187) 


Note that the standardized gaussian random vector f, is the same at the predictor and corrector 
levels. This is to ensure that the numerical approximation given by Eqs. (185)-(187) reflects the 
Markovian behavior of the general diffusion process. 96 Markovian or “memoryless” processes are 
stochastic processes in which future states are not influenced by past behavior. 87,97,89 

The drift coefficients D, and B require knowledge of the mean field velocity and viscosity. Inter- 
polation is required for these quantities since the Lagrangian particles are not restricted to the 
finite difference grid points. Fourth order Lagrange polynomials are used to interpolate the desired 
quantities from the grid to the particle location. 

Each particle contains information regarding the composition of the scalar field. This includes 
the density weighted species mass fractions (and temperature if the chemistry model requires it). 
Let denote the value of the a t>l scalar (o = 1,2, .... ) for the k ih particle ( k — 1,2,..., A) 

located at the Lagrangian coordinate X*’ as described in section 3.3.4. At each time step these 
compositional values are subject to change due to the effects of molecular mixing and chemical 
reaction. Using Dopazo’s relax to mean model for molecular mixing the composition of each 
particle changes according to 

(**)”"'* = (4>i) n - ((</>o) n - <<^>)exp[-icV t A t] + <4fe (188) 

at each timestep. The turbulent mixing frequency u t is backed out from the turbulent viscosity 
and length scale using the relation w, = 3;y/A 2 . After the mixing step has been completed the 
particles undergo reaction. This is performed by sweeping over all the particles and determining 
the fine grain reaction rates u For example, for the simple reaction A + B — ♦ V with reaction 
rate u>a = -kj{pfA)(p/B)> the fine grain reaction rates are given by 

* k A = = -W^aKpH)- ( 189 ) 
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The new particle composition is then determined from: 


(4> k a ) n+1 = (<P k a ) m,x + u h Q At. 


(190) 


3.4.3 Determination of Filtered Quantities 


For homogeneous flows averages are determined by simply summing over all the particles: 


N N 

<A(4>)> = J2M4> k ) = J2 Ak - 


k=l k = 1 


(191) 


For inhomogeneous flows the situation is complicated by the fact that the PDF varies spatially. A 
discrete summation consistent with the Large Eddy PDF is: 


-r, _ EfcLi>l(X*(0,0G(X*(*)-x) _ Y.k=i A^jx) 

Ete! G(X*(<)-x) ELi^-(x) 


(192) 


Averages and higher order moments can be calculated in this manner. In the present work where 
a uniform filter of width A = 2Ax = *2A y has been used, calculation of averages at the grid point 
with coordinates (i,j) reduces to summing over all particles in the square region of dimension 2Ax 
by 2 A y centered at the grid point. 


For Reynolds averaged solutions, similar procedures must be taken to generate meaningful ensem- 
bles. However it should be noted that in the case of LES the spatial dependency of the ensemble 
is reflected by the filter G(x) in the definition of the large eddy PDF (see Eq. (3.3.5)); for the case 
of RANS there is no such spatial dependency in the “conventional PDF”. Rather, the practice of 
constructing an ensemble out of particles within some volume of a point is out of necessity since in 
general the Lagrangian particles will not coincide with the Eulerian grid points. In this sense LES 
is more in tune with the Lagrangian solution procedure. 


3.5 Results 

To demonstrate the feasibility of LES-PDF, one non-reacting and one reacting simulation w f ere 
performed. For comparison, two additional runs were performed in which the Favre averaged 
species equations were solved by the finite difference method described in section 3.4.1 rather than 
the LES-PDF methodology. For the reacting finite difference simulation the mean reaction rate 
was modeled as uJ a = uj a (<j>). It is well recognized that such an assumption may be ofT by orders 
of magnitude and gross errors can be incurred. This simulation was performed to compare to the 
LES-PDF reacting simulation in which the effects of reaction are accounted for exactly. In order 
to provide meaningful results, an additional control simulation was run using the LES-PDF using 
= w Q (<£) to elucidate any discrepancies with the finite difference procedure. 
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All simulations were performed at a Reynolds number of 800 based on initial vorticity thickness 
Su/O’ The physical domain in all cases is 60^0 by 306^ which was discretized into a 256 by 128 
computational domain. The velocity profile at the inlet consists of a 4 to 1 velocity ratio with a 
hyperbolic tangent distribution. Forcing of the cross stream velocity at the frequency corresponding 
to the most unstable mode for the hyperbolic tangent velocity profile was used to perturb the layer 
and induce coherent large scale structures. 98 Zero first derivative conditions are assumed at the 
freestreams. The chemistry model utilized is an irreversible second order isothermal reaction of the 
form A + B — ► V. In each of the reacting simulations the Damkohler number Da = kj /[U r€ j /6 U0 \ 
was set to 2. 

In order to save time during the computation, only the region surrounding the reaction zone from 
y = —3.75 <5 Wo to y = 3.75 b ^ 0 is initialized with the Monte Carlo particles (one quarter the com- 
putational domain). In each cell in this region 30 particles were randomly placed. Thus on average 
roughly 120 particles fall within the filter width 2Ar by 2Ay to constitute an ensemble, although 
this number may vary somewhat as particles convect and diffuse from cell to cell. Additionally 
as particles exit the domain, new particles are introduced at x = 0 with a randomly chosen y 
coordinate in the vicinity of the reaction zone. 

Figure 3.1 is a contour plot displaying the particle density. Note how the Monte Carlo particles 
are only distributed in the region surrounding the reaction zone. There is some variation locally, 
however overall the particle density is roughly uniform at a value of 120 particles per ensemble. 
The oscillatory behavior of the particle zone in the last third of the domain is suggestive of some 
organized large scale structure. Indeed such coherent structures are evident in Fig./ 3.2 which 
compares contours of species A Favre averaged mass fraction for LES-PDF (3.2(a)) to those of 
standard finite difference LES (2 — b) for the non-reacting case. The solution provided by the PDF 
methodology compares favorably with the standard procedure. Such comparison is important in 
that it demonstrates that the particle method is capable of resolving the convective and diffusive 
transport mechanisms. 

Figure 3.3 displays the contours of species A mass fraction for the reacting case. Although at first 
glance the plots appear similar, it is soon apparent that the gradients in Fig. 3.3(b) (finite difference 
with & Q = u Q (<}>)) are steeper suggesting the reactants are separated by a relatively thin reaction 
zone. This is indicative that the reaction is more vigorous and has progressed further compared to 
the LES-PDF solution depicted in Fig. 3.3(a). This is more clearly seen in the product mass fraction 
contours in Fig. 4. Substantia] less product formation is predicted by the LES-PDF simulation and 
the solution is significantly more diffuse. Clearly this is a result of the LES-PDF solution resolving 
the SGS subgrid fluctuations while these terms are neglected in the finite difference solution. 

To assist in further appreciation of the effects due to the SGS scalar fluctuations the covariance 
(p/aY^/bY ^ plotted in Fig. 3.5(a) for the non-reacting case and Fig. 3.5(b) for the reacting 
case as determined by the solution to the large eddy PDF. The negative values of this quantity 
as predicted by the LES-PDF procedure indicate that neglect of this subgrid term would result in 
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unphysically high reaction rates and product conversion. Further evidence is given by reaction rate 
contours predicted by LES-PDF in Fig. 3.6(a) and that predicted by the finite difference solution 
in Fig. 3.6(b). Clearly the reaction rates in the LES-PDF solution are an order of magnitude 
lower than that predicted by the solution neglecting the SGS scalar fluctuations. This indicates 
the approximation Zo a = io a (4>) cannot be justified. While the scalar covariance is very difficult to 
model, LES-PDF avoids this obstacle altogether as it has the distinction of being able to resolve 
this term. 

To be confident that the LES-PDF calculations are representative of the true physics, an additional 
simulation was performed with the particle method in which the reaction rate was improperly 
“modeled” by ZJ Q — u> a {4>)- Figure 3.7 is a contour plot of the product V mass fraction. Notice 
that the product distribution exhibits a remarkable likeness to the finite difference solution to the 
species transport equations with the scalar SGS terms neglected. Whether the PDF method or 
the finite difference procedure is used the product formation is grossly overpredicted if the subgrid 
correlations are neglected. Figure 3.8 displays the reaction rate contours for this simulation. As 
expected the predicted reaction rate is an order of magnit ude higher due to the neglect of the subgrid 
terms. This agreement between the finite difference and the PDF procedures is important since it 
elucidates that the neglect of the SGS terms rather that the difference in numerical procedure is 
responsible for the discrepancy of the previous runs. 


3.6 Conclusion 

A PDF method suitable for chemically reactive flows is developed in the context of large eddy 
simulation. The clear advantage of PDF methods is their inherent ability to resolve SGS correlations 
that otherwise have to be modeled. Because of the lack of robust models to accurately predict these 
correlations in turbulent reactive flows, simulations involving turbulent combustion are often met 
with a degree of skepticism. The PDF methodology avoids the closure problem associated with 
these terms but rather treats the reaction exactly. 

The first LES-PDF simulations of a chemically reactive flow have demonstrated the feasibility of 
utilizing Monte Carlo methods in the context of LES. Comparison with a finite difference solution 
which does not attempt to model the SGS scalar covariance indicates that neglect of this term leads 
to unphysically high product conversion rates. 

While the present work indicates there is much promise in the LES-PDF methodology, much work 
needs to be done. Utilization of compressible dynamic SGS models is a natural starting point. 
Furthermore, the development of LES with variable grid and filter spacing as addressed by" is an 
important step towards the ability of such models to predict engineering flows with complex geome- 
tries. PDF transport methods require the closure of molecular mixing; in the present research this 
has been addressed with a simple deterministic model. Mixing models require a turbulent mixing 
frequency. Proper determination of this mixing frequency is desirable. Additionally, inclusion of 
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the energy in the compositional domain is a necessary step in order to simulate reactive flows with 
temperature dependent kinetics. Furthermore, since the reaction is represented without approxi- 
mation realistic chemical kinetics present no difficulty to the particle method and hence should be 
investigated. Much of this work is already underway. 

Additional extension of the LES-PDF methodology to solve the joint velocity-composition PDF is 
also an area to be explored. This area has seen considerable development in the area of RANS by 
Pope. 100 The velocity-composition PDF has the added advantage that velocity correlations appear 
in closed form and eddy viscosity models are not required to model their effects. 72 In the case 
of LES, however, it is expected that the subgrid models are more universal in their applicability 
and hence the dependency on SGS closures for the scalar PDF method does not pose too much 
concern. Furthermore, particle methods as described in this report can be incorporated into the 
large stockpile of existing finite difference and finite volume CFD codes. 


4 Work in Progress 


We are currently in the process of combining the techniques discussed in the last tw'o sections for 
LES of high speed reacting turbulent flows. The primary tool in this part of our activities is the 
Monte Carlo LES-PDF solver as discussed in the last section. However, our aim is to replace 
the Smagorinsky hydrodynamic subgrid model with an algebraic closure. For this purpose, the 
same computer code used for LES-PDF is being modified, presently, we are only considering LES 
of incompressible reacting flows. Therefore, only the models developed previously 19 ' 1 are being 
utilized (of course with proper modifications to make them suitable for LES). At this point, there 
are some computational problems that need to be resolved. We hope to have some results before 
our next report is due. 

In the upcoming year (Year 3 of phase II), our efforts will be concentrated on the following tasks: 

(1) Completion of the mathematical formulation and computational implementation of LES via 
algebraic closures. 

(2) Generation of extensive computational results via RANS of high speed flows with and without 
chemical reactions and comparisons with experimental data. 

(3) Fine tuning of the LES-PDF code. Application of this code for simulations of several reacting 
flow configurations. Comparisons with experimental data. 

(4) Extension of LES with consideration of the joint PDFs of subgrid velocity and scalars. Our 
efforts will be concentrated primarily on the mathematical formulations and related issues. 
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5 Personnel and Management 


Dr. Peyman Givi is the PI of this project. Drs. Dale B. Taulbee and Cyrus K. Madnia serve as the 
Co-PIs. As before, Dr. Givi is responsible for a timely and successful completion of all the tasks. 

There are presently two Ph.D. candidates who are being supported full-time by this grant: Mr. 
Virgil Adumitroaie and Mr. Paul J. Colucci. Mr. Adumitroaie has been supported under this 
program since the initiation of Phase II (and departure of Dr. Steven H. Frankel). Mr. Colucci is 
replacing a former RA, Mr. Craig Steinberger. The reason for this replacement is that the work of 
Mr. Colucci is more closely related to this program. Mr. Steinberger is continuing his work in DNS 
of reacting flows the results of which will be utilized for this work. However, he is no loner being 
financially supported under this grant. 
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Figure Captions 


Figure 2.1: The reduction of the averaged residual for the Gottlieb-Turkel scheme. 

Figure 2.2: Downstream evolution of the shear-layer width at steady-state. 

Figure 2.3: Cross-stream variation of ( U - U\)/(U\ - U 2) 2 for the mixing layer. 

Figure 2.4: Cross-stream variation of {u 2 )f(U\ — U 2) 2 for the mixing layer. 

Figure 2.5: Cross-stream variation of (u'v f )/(Ui - U 2) 2 for the mixing layer. 

Figure 3.1: Particle number density (per ensemble) contours . 

Figure 3.2: Contour plots of species A mass fraction (Da=0): (a) Monte-Carlo LES-PDF, (b) finite 
difference. 

Figure 3.3: Contour plots of species A mass fraction (Da=2): (a) Monte-Carlo LES-PDF, (b) finite 
difference. 

Figure 3.4: Contour plots of product V mass fraction (Da=2): (a) Monte-Carlo LES-PDF, (b) 
finite difference. 

Figure 3.5: Contour plots of SGS species covariance : (a) Da=0, (b) Da=2. 

Figure 3.6: Reaction rate contours: (a) Monte-Carlo LES-PDF, (b) finite difference. 

Figure 3.7: Product V mass fraction contours. Results generated by Monte Carlo LES-PDF 
assuming £>(<£) = cj(0). 

Figure 3.8: Reaction rate contours. Results generated by Monte Carlo LES-PDF 
assuming U (</>) = Cj{4>), 
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Figure 1: Particle number density (per ensemble) contours 
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Figure 3.3: Contour plots of species A mass fraction 
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Figure 3.4: Contour plots of product P mass fraction 

(Da.=2) (a) Monte-Carlo LES-PDF, (b) finite difference 
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Figure 3.6: Reaction rate contours: (a) Monte-Carlo LES-PDF, 
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Figure 3.7: Product P mass fraction contours Results 

generated by Monte -Carlo LES-PDF assuming <AB>=<A><B> 



Figure 3.8: Reaction rate contours. Results generated 

by Monte -Carlo LES-PDF assuming <AB>=<AxB> 








